arXiv: 1503.00023vl [astro-ph.SR] 27 Feb 2015 


Mon. Not. R. Astron. Soc. 000 , 000-000 (0000) Printed 3 March 2015 (MN IAT^X style file v2.2) 


The MLP Distribution: A Modified Lognormal Power-Law 
Model for the Stellar Initial Mass Function 


Shantanu Basu 1,2 *, M. Gil 2 f and Sayantan Auddy 1 

1 Department of Physics and Astronomy, The University of Western Ontario, London, ON, Canada N6A 3K7. 

2 Department of Applied Mathematics, The University of Western Ontario, London, ON, Canada N6A 5B7. 


3 March 2015 


ABSTRACT 

This work explores the mathematical properties of a distribution introduced by Basu 
& Jones (2004), and applies it to model the stellar initial mass function (IMF). The 
distribution arises simply from an initial lognormal distribution, requiring that each 
object in it subsequently undergoes exponential growth but with an exponential dis¬ 
tribution of growth lifetimes. This leads to a modified lognormal with a power-law 
tail (MLP) distribution, which can in fact be applied to a wide range of fields where 
distributions are observed to have a lognormal-like body and a power-law tail. We 
derive important properties of the MLP distribution, like the cumulative distribution, 
the mean, variance, arbitrary raw moments, and a random number generator. These 
analytic properties of the distribution can be used to facilitate application to modeling 
the IMF. We demonstrate how the MLP function provides an excellent fit to the IMF 
compiled by Chabrier (2005) and how this fit can be used to quickly identify quantities 
like the mean, median, and mode, as well as number and mass fractions in different 
mass intervals. 
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1 INTRODUCTION 

The distribution of stellar and substellar masses at birth, 
the Initial Mass Function (IMF), is a key feature of star for¬ 
mation. It has been studied intensely since first estimated 
by Salpeter (1955), who measured a power-law tail for high 
masses of the approximate form dN/d\n M oc M -1 ' 35 . Sub¬ 
sequent work has established a shallower slope at masses 
less than and a turnover at approximately 0.2 Mg when 
the masses are placed in logarithmically spaced bins. These 
compilations of the IMF have tended to identify a power- 
law profile at high masses (e.g., Scalo 1998; Kroupa 2001, 
2002), although earlier work (Miller & Scalo 1979) did fit 
the IMF with a lognormal distribution. Chabrier (2003) (see 
also Chabrier (2005)) has compiled an IMF in the sub- 
stellar and low mass stellar regime and advocates a log¬ 
normal fit for masses up to Mq and a power-law fit for 
M > Mq. By appealing to the Central Limit Theorem 
(CLT), Chabrier (2003) claims a better rationale for the 
lognormal fit at low masses than the approach of using bro¬ 
ken three-component power-laws (Scalo 1998; Kroupa 2001, 
2002). However, Chabrier’s overall fit also requires joining 
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a lognormal with a power law at high masses, so has one 
joining condition instead of two. An ideal next step is to 
find a single function with no joining conditions that has a 
rationale that is at least on par with appeals to the CLT. 

There are many disciplines in which a desired distribu¬ 
tion is one that is like a lognormal at low and intermediate 
values, with a characteristic peak and turnover, but tran¬ 
sitions to a power law distribution at high values. Besides 
astronomy, this need has arisen in fields as diverse as biol¬ 
ogy (Limpert et al. 2001), computer science (Mitzenmacher 
2006), ecology (Allen et al. 2001), and finance (Mandel¬ 
brot 1997), to name several. A review of common statis¬ 
tical resources reveals that very few analytic functions of 
this type exist, hence the attempts to fit empirical distribu¬ 
tions by patching together different functions over different 
domains. Power-law distributions were first introduced by 
Pareto (1896) to explain the distribution of incomes seen 
in data from many different countries. The distribution of 
city sizes also shows a power-law character and regularity 
across many countries. It is referred to as Zipf’s Law (Zipf 
1949). Henceforth, we refer to pure power-law distributions 
as Pareto distributions, in conformity with much of the sta¬ 
tistical literature. 

In this paper, we analyze and characterize the proper¬ 
ties of a hybrid three-parameter probability density function 
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(pdf) introduced by Basu & Jones (2004). We feel that it 
can be used to fruitfully model data sets that exhibit both 
lognormal-like and power-law behavior. Indeed, the MLP 
distribution illustrates the fact that many generative pro¬ 
cesses that lead to a lognormal distribution, can with some 
modification, yield a power-law tail instead. As a result of its 
origin as a modified lognormal, two parameters of the MLP 
distribution are identified as po and ag, preserving the nota¬ 
tion of the lognormal distribution, while the third parameter 
is a, the power-law index that also characterizes the Pareto 
distribution. However, the parameters fig and ag no longer 
represent the mean and variance, respectively, of the loga¬ 
rithm of the random variable, as they do for the lognormal 
distribution. This three-parameter function is a simpler and 
more readily usable form of a more general four-parameter 
pdf derived by Reed in several papers (Reed 2002; Reed 
2003). The latter function arises from a stochastic growth 
law, Geometric Brownian Motion, rather than the pure ex¬ 
ponential growth used by Basu & Jones (2004). An advan¬ 
tage of the pdf we use here is that it introduces only one 
additional parameter beyond that in the lognormal. As a 
result, it is a natural first step when fitting data that may 
look like a modified lognormal, and its relatively compact 
analytic closed-form expression makes it easy to use with 
common fitting techniques. 

The model of Basu & Jones (2004) falls into one of two 
major categories of IMF models in modern day star forma¬ 
tion theory, that of an accretion based scenario in which 
temporal effects are important. The basic idea is that star 
formation is a killed (or stopped) process, with mass accre¬ 
tion terminated by events such as stellar outflows, ejection 
from the mass reservoir, or any other reason for emptying the 
mass reservoir. The distribution of accretion lifetimes then 
plays a key role in setting the shape of the IMF, while the 
Jeans mass does not. Other models in this category include 
those of Adams & Fatuzzo (1996), Bate & Bonnell (2005), 
and Myers (2009, 2014), with the latter developing a model 
that can fit the observed IMF by tuning several parame¬ 
ters. The alternate scenario is one in which the mass dis¬ 
tribution is determined by the spatial properties of the gas 
distribution and by gravitational instability, so that some 
combination of the turbulent spectrum and the Jeans mass 
set the IMF, e.g., Padoan & Nordlund (2002); Hennebelle 
& Chabrier (2009), with the latter developing an analytic 
approach. Although the fitting of analytic functions to the 
IMF cannot alone settle which models may be most suit¬ 
able, they can however greatly facilitate analysis in fields 
like galaxy studies, where the conclusions depend strongly 
on an adopted IMF model. Analytic functions also allow a 
simpler analysis of the effect of varying IMF parameters. 

The paper is organized as follows. In Section 2 and Sec¬ 
tion 3 we introduce the lognormal and Pareto distributions, 
respectively, and present some of their relevant properties. 
This is done for completeness of the presentation and to add 
context when reading the new results about the MLP dis¬ 
tribution. In Section 4 we discuss the formulation of Basu & 
Jones (2004) that leads to the MLP distribution. We then 
examine some relevant mathematical properties of this dis¬ 
tribution in Section 5, which includes expressions for its cu¬ 
mulative distribution function, mean and variance, arbitrary 
raw moments, and an approximation to its mode. These ex¬ 
pressions, excluding the approximate mode, are shown to 


reduce to the corresponding lognormal expressions in the 
appropriate limit. In Section 6, we fit the MLP distribution 
to the IMF of Chabrier (2005) and then use the analytic 
function to quickly estimate some of the above described 
IMF properties, as well as a cumulative mass fraction. Some 
closing remarks are given in Section 7. The derivations of the 
expressions given in Section 5 are included in Appendix B, 
which also contains relevant properties of the error and com¬ 
plementary error functions and related integrals used in this 
paper. 


2 THE LOGNORMAL DISTRIBUTION 

According to the CLT of probability theory (Gut 2005), 
if Xi,...,X n are identically distributed, independent ran¬ 
dom variables with mean p, and standard deviation a, then 
Z = Q ZXi — nn) /(y/na) converges in distribution to the 
standard normal variable N( 0,1) with zero mean and vari¬ 
ance of unity. In addition, the identical distribution assump¬ 
tion for Xi, ...,X n may be dropped, and the result will also 
follow provided certain conditions (Lindeberg’s) are satisfied 
(Gut 2005). As pointed out in Golberg (1984) the CLT is 
used to partially justify why so many observable phenomena 
appear to be normally distributed: if the variable of interest 
is thought to be influenced by the sum of a large number 
of independent factors, then the CLT can be invoked to ex¬ 
plain the apparent normality. However, it is frequently the 
case that the variables of interest are non-negative, whereas 
normal variables are not. One way to circumvent this com¬ 
plication is to have a distribution with similar properties 
but which is always positive. Suppose that Z above was in¬ 
stead assumed to be of the form Z = Xi ■ Xn ■ ... ■ X n . Then 
W = In Z = where Y, = lnX*. Note that Y u ...,Y n 

satisfy the conditions of the CLT, so that W converges to 
a normal random variable, W ~ N(n,a 2 ). Moreover, if W 
is normal, then Z = e H has a lognormal density (Golberg 
1984; Aitchison & Brown 1957; Crow & Shimizu 1988), 

Thus, under the assumptions above, Z converges to a lognor¬ 
mal. We list some of the relevant properties of the lognormal 
distribution below: 


(i) Raw Moments: 


E[Z \ = exp ( k/j, + -k ct 


(ii) Mode: 

z m = exp(n — a 2 ) ; 

(iii) Variance: 

Var(Z) = exp (2 fi + a 2 ) [exp(<r 2 ) — 1] ; 

(iv) Cumulative Distribution Function (cdf) 

'\nz — 


F z (z;fi,a ) = - 


1 + erf 


V2o 


( 2 ) 

( 3 ) 

( 4 ) 

( 5 ) 


where erf(z) is the error function (see Equation Al in the 
appendix). 
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3 THE PARETO DISTRIBUTION 

Pareto distributions have been used extensively to model a 
wide variety of phenomena in the sciences and social sci¬ 
ences, such as the size of forest fires, the intensity of earth¬ 
quakes, the citations to papers, and the population of cities. 
For a recent review, see (Clauset et al. 2009). 

A random variable X has a Pareto distribution if its 
probability density function is 


/(*) = ’ * > b > 0 ( 6 ) 

(Pareto 1896; Mitzenmacher 2006; Seggern 1990), where a > 
0. We list some of its more relevant properties below: 

(i) Raw Moments: 

= q, > „ ; ( 7 ) 

a — n 

(ii) Mode: 

Xm — b , ( 8 ) 


(iii) Variance: 

2 ab 2 

a (a — l) 2 (a — 2) ’ 2 ’ 

(iv) Cumulative Distribution Function (cdf) 


(9) 


F x (x; a,b) = 1 — — 

. x 


( 10 ) 


4 THE MLP DISTRIBUTION 


The key element in the derivation of the MLP distribution is 
that even though the initial values of a random variable may 
have a lognormal density, later time evolution could in fact 
skew their distribution. In other words, the subsequent time 
of growth of the random variable (representing a physical 
quantity) is itself another random variable. Our pdf can be 
derived analytically on this basis using a few simplifying 
assumptions. 

Imagine that initial values of a quantity Mo are drawn 
randomly from a lognormal distribution with parameters no 
and tJo. If the subsequent growth of an object with Mo = mo 
is characterized by exponential growth 


dm 

dt 


7 m, 


( 11 ) 


with fixed growth rate 7 , then the multiplicative relation 
between values m of individual objects is preserved. A log¬ 
normal distribution is then maintained with the same ao but 
the mean of the logarithmic values shifts to n 0 + yt after a 
fixed time t. However, we can treat the time as a random 
variable and draw it from an exponential pdf 

f(t) = Se~ st , ( 12 ) 


where 5 is a stopping rate. In this case, we can derive the 
probability density function of hnal masses m as 



8e~ st 

y/Tnaom 


exp 


[lnm — no ~ 7 1] 2 
2 °'o 


dt . (13) 


Using the identity in Equation A9 we find the closed form 
f(m) exp (a/ito + a 2 al/2) m~ (1+a) 

, f 1 f In m- no \\ r r. , 

x erfc — 1 = aao -, m £ 0, 00 ) , 

\ v 2 V °o )) 

(14) 

where a = 5/y. 

Figure 1 shows the MLP density /(m) (Equation 14) 
and a lognormal distribution for similar values of param¬ 
eters. The specific values of the parameters used here are 
somewhat arbitrary, but correspond approximately to a best 
fit lognormal for the low mass end of a stellar mass distri¬ 
bution, as modeled in Basu & Jones (2004). While no is 
largely a scale-dependent parameter, the other two parame¬ 
ters are expected to fall in the approximate range 0 ^ 07 ^ 1 
and 1 ^ a ^ 2 based on fits of lognormal and power- 
law distributions to a wide range of phenomena in the 
sciences and social sciences (Limpert et al. 2001; Clauset 
et al. 2009). The function derived here is related to a four- 
parameter distribution (Reed 2002; Reed 2003) that can be 
derived under the assumption of geometric Brownian mo¬ 
tion, dm = ymdt+amdw, which is a stochastic growth law 
where dw represents J?(0,1), the uniform random variate in 
the interval [0, 1], and a is an amplitude of the fluctuations. 
The resulting pdf is a double-tailed Pareto distribution, with 
coefficients that are roots of a quadratic equation. The pdf 
we derive here results from a simpler model and has the ad¬ 
vantage of a single-expression closed form. It may also more 
easily correspond to many data sets that seemingly warrant 
only one additional parameter (beyond the original two of 
the lognormal distribution) in order to get a reasonable fit 
and quantify the data. 


5 RELEVANT MATHEMATICAL 

PROPERTIES OF THE MLP DISTRIBUTION 

A few relevant mathematical properties of the MLP distri¬ 
bution are examined below, leaving most of the necessary 
calculations in Appendix B. Let M denote a random vari¬ 
able with an MLP distribution of density f(m). 


5.1 Cumulative Distribution 


The MLP cumulative distribution function, 

/ m 

f(t)dt , (15) 

- OO 


is given by 
F(m) — ^erfc ( — 


In(m) - no 
\[2oo 


- - exp ( ano + 


2 2 
a a 0 


"erfc 


/ aero 


ln(m) — no \ 
p2ao ) 
(16) 


Notice that as m —> 0 we can see by the result in Equa¬ 
tion B3 that the behaviour of F(m) is dominated by the first 
term, which is exactly the cumulative distribution function 
for a lognormal random variable (Equation 5) with param¬ 
eters no an d <r 0 . 
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m 



Figure 1. Comparison of a lognormal density function with [i = 
— 1.5 and a = 0.5 with the MLP density function with /xo = —1.5, 
ctq = 0.5, and a = 1. The mass is normalized by a solar mass Mq. 


5.2 Mean, Variance, and Raw Moments 

The Arth raw moment of M, defined as the expectation value 
of M k , 

r oo 

E[M k ] = / m k f(m)dm, (17) 

■J o 

exists if and only if a > k, in which case it is given by 

E[M k ] = a exp f - + fj, 0 k) , a> fc. (18) 

Note that this expression is exactly the formula for the raw 
moments of a lognormal distribution (Equation 2) with pa¬ 
rameters fio and fro, scaled by the factor a/(a — fc), and in 
the limit as a —> oo, the expressions are identical. This is 
consistent with the derivation of the MLP distribution in 
Section 4. The limit a —> oo corresponds to 7 —> 0 for finite 
death rate 5, so that the drift term vanishes, and the distri¬ 
bution remains a lognormal with mean fi = hq. Assuming 
a > fc, we can obtain the following expressions, including 
the mean and variance of the distribution: 

E[M] = -^-j- exp + /to j , a > 1, (19) 


E[M 2 ] = exp (2 [al + H) , « > 2, (20) 


Var(M) = E[M 2 ] - (E[M]) 2 


= a exp(oo + 2 /xq) 


:-2 {a-If 


a >2. 


( 21 ) 


Higher moments around the mean can be computed using 
Equation 18 with the identity 


3 =0 


E[(M - E[M]) n ] = ^2 ( ) E[M j ](-l) n ~ j {E[M]f 


( 22 ) 


5.3 Mode 

To find the mode, that is, the value m* that maximizes the 
MLP pdf in Equation 14, we must solve the transcendental 
equation 

f'(m) = 0 4=4> Aerfc(-u) = e~ u , (23) 

where 

y r , fn If lnm—no\ 

K = cro(a + 1 ) 1 / 77 , u=—=latro - . (24) 

V ^ v2 \ cro / 

Although the solution to Equation 23 will generally require 
the implementation of numerical methods, we note that if 
Km 1 then u = 0 provides an approximate solution to 
Equation 23, which in terms of the original parameters re¬ 
sults in 

m* = exp (/r 0 + aoo) . (25) 

The approximation in Equation 25 is useful only when the 
assumption K m 1 is closely met, and behaves poorly oth¬ 
erwise. However, when a precise numerical solution is re¬ 
quired, one can use this approximation as a starting point 
in the iteration procedure being implemented. It is worth 
noting that even if one tried to find the peak in the space 
mf(m) vs log(m), the resulting equation to be solved is in 
fact the same, owing to the fact that mf(m) is just f{m) 
with a different value of a, and d/d(\og(x)) = x(d/dx). 


5.4 Random Number Generation 

For practical purposes of comparing data sets with a model 
distribution, it is valuable to be able to draw a random sam¬ 
ple from the model. Using the definition of the lognormal 
random variable (see Section 2), a random number drawn 
from a lognormal pdf (Equation 1) is an element of the log¬ 
normal random variate 

L(/i,a) ~ exp(/r + aN(0, 1)), (26) 

where 1V(0,1) is the normal random variate with zero mean 
and variance of unity. For drawing from an exponential pdf 
(Equation 12), the exponential random variate is 

E{6) ~ — .T 1 ln(R(0,1)), (27) 

where i?(0,1) is the uniform random variate in the interval 
[0,1]. The above formula can be obtained from Equation 12 
through the general method of calculating the cdf, inverting 
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it, and then drawing the argument as an element of the uni¬ 
form random variate R(0, 1). We can use Equation 26 and 
Equation 27 to derive a random variate for the MLP distri¬ 
bution. We note that the MLP distribution (Equation 14) is 
formally obtained from an initial lognormal pdf with mean 
Ho under the transformation ho —> Ho + "ft, in which t is 
chosen randomly from an exponential distribution. There¬ 
fore, we can use Equation 26 and Equation 27 to write that 
the MLP random variate is 

M(ho, ct o, a) ~ exp(Ho+&o N( 0,1) — a -1 ln(i?.(0,1))). (28) 

Equation 28 can be used to draw a random sample from 
the MLP pdf just as Equation 26 can be used to draw a 
random sample from a lognormal pdf. 


6 FITTING THE MLP FUNCTION TO THE 

IMF 

Here, we find a set of parameters for the MLP distribu¬ 
tion that fit the updated IMF of Chabrier (2005). We use 
a normalized form of Equation (1) in Chabrier (2005). It 
is divided into logarithmic bins of mass normalized by Mg 
with intervals of 0.05 over the range of 0.06 — 100 M 0 and 
we use a Levenberg-Marquardt least-squares minimization 
method to fit the MLP function. The best fit parameters 
are Ho = —2.404, <jo = 1.044, and a = 1.396. Figure 2 shows 
the best fit MLP function along with the Chabrier func¬ 
tion, both in normalized form. Figure 3 shows three sets of 
random samples drawn from this best fit MLP function, us¬ 
ing Equation 28 and sample sizes of 100, 1000, and 10000, 
respectively. The samples are made into histograms with 
logarithmic bin width of 0.2, and each histogram is repre¬ 
sented by a different symbol. We find that samples of well 
over 100 are needed to adequately populate the power-law 
tail. For MLP samples of about 100 or less, fitting a log¬ 
normal pdf typically yields a better goodness-of-fit statistic 
than the MLP distribution, given that it also has one less 
fitting parameter. 

Some sample routines for fitting the MLP function 
or drawing random samples from it can be found at 
http: //www. astro .uwo. ca/~basu/mlp .html. 

With a best fit MLP function in hand, it is relatively 
easy to quantify several important properties of the IMF. 
Equation 19 can be used to calculate the mean mass in the 
distribution; it is 0.55 M 0 . If we were to truncate the dis¬ 
tribution at 100 M 0 , the mean mass is about 10% lower 
at 0.49 Mg. The mode of the distribution f(m) is calcu¬ 
lated numerically to be 0.04 M 0 . However, the usual prac¬ 
tice is to plot data as dN/dln M = mf(m) rather than 
dN/dM = f(m), so the usually quoted IMF peak is that of 
mf(m .). We find this to be 0.16 Mg for our best fit param¬ 
eters. The mode of mf(m) is found analogously to that for 
f(m) described in subsection 5.3, but with a lowered by one 
integer value. 

Figure 4 shows the cumulative MLP distribution for the 
best fit parameters. The legend panel in Figure 4 illustrates 
the specific values for F(m) for several relevant mass values, 
e.g., one solar mass, as well as the minimum stellar mass, 
m = 0.075M 0 (Chabrier & Baraffe 2000). Although the data 
used to generate the IMF quoted by Chabrier (2005) does 
not span the entire hypothetical mass range m £ [0, oo], 



log(m) 

Figure 2. Comparison of the Chabrier IMF with its correspond¬ 
ing best fit MLP, where ho = —2.404, rro = 1.044, and a = 1.396. 
The mass is normalized by Mg. 


one can use the numbers in Figure 4 to quickly estimate 
some highlights of the mass function, with interesting pairs 
of m, F(m) marked by symbols in Figure 4: about one quar¬ 
ter of objects are substellar, about one third of objects have 
masses less than 0.1 Mg, the median mass is 0.17Mg, just 
over 90% of objects have masses less than 1 Mg , and only 
9% of objects have masses in the range 1 — 10Mg. Fig¬ 
ure 5 shows the mass fraction 4>(m) = m' f(m') dm' 
(with closed form given in Equation B15) for the best fit 
parameters, normalized by the value at m = 100Mg. Sym¬ 
bols denote interesting pairs of m, c/>(m), and many provide 
an interesting contrast to their counterparts in Figure 4. For 
example, only about 2% of the total mass is in substellar ob¬ 
jects, a majority of mass is tied up in objects more massive 
than lMg, and half of all mass is in objects above mass 
1.4 Mg. 

The ability to model stellar populations with an IMF 
that has clearly identifiable number or mass fractions in dif¬ 
ferent mass ranges should be a key aid to galaxy studies 
(e.g., Kennicutt 1998), for example in the determination of 
star formation rates, chemical evolution, and mass-to-light 
ratios. 


7 CONCLUSIONS 

We have derived several important properties of the mod¬ 
ified lognormal power-law (MLP) probability distribution 
function that has been recently introduced in the literature. 
The three-parameter MLP function has the salutary prop¬ 
erties of a main body that resembles a lognormal, including 
a peak value and a decline toward low values, as well as a 
power-law tail at high values. This function can potentially 
be applied in a variety of fields where empirical distribu¬ 
tions may have a power-law tail that coexists with a peak at 
lower values. The MLP distribution can also help to settle a 
frequent contentious question: is a data set consistent with 
a lognormal distribution or does it also show evidence for 
a power-law tail? Simultaneous fitting of the lognormal and 
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Figure 3. The MLP function wth /xo = —2.404, cr 0 = 1.044, and 
a = 1.396, overlaid with histogram values for random samples 
drawn from the distribution with size N = 100 (blue circles), 
N = 1000 (black squares), and N = 10000 (green triangles). 
All histograms are binned in increments Alogm = 0.2, and his¬ 
togram values are the fractional number in each bin divided by 
A lnm. 



log(m) 

Figure 4. MLP cumulative distribution with //q = —2.404, a o = 
1.044, and a = 1.396. The mass is normalized by Mq. The legend 
illustrates interesting values of the pairs (m, F(m)). 

MLP distributions to the same data sets can help to answer 
this question in many fields. 

Comparison of real data sets with the MLP distribution 
is facilitated by the results presented in this paper. We have 
derived analytic expressions for the cumulative distribution 
function, the mean, variance, and higher moments of the dis¬ 
tribution. We also derived an approximation for the mode 
that can serve as an initial guess for nonlinear solvers that 
can iterate to a more exact solution. The random variate of 
the MLP distribution has also been introduced. Together, 
these results can put the MLP distribution on a more equal 
footing with many classical distributions (e.g., lognormal, 
exponential, Pareto, Rayleigh) that are frequently used to 
fit empirical data. The use of the MLP distribution to fit 


log(m) 

Figure 5. Mass fraction for an MLP distribution with //.q = 
—2.404, (To = 1.044, and a = 1.396. The mass is normalized by 
Mq. The legend illustrates interesting values of pairs of m and 
the mass fraction. 

an empirical data set has been demonstrated and we show 
how useful information about the mean, mode, median, and 
distributions of number or mass fractions can be easily cal¬ 
culated. 


8 ACKNOWLEDGEMENTS 

SB acknowledges the hospitality of the Isaac Newton In¬ 
stitute of Mathematical Sciences at Cambridge University 
during the initial stages of writing of this paper. This work 
was supported by an NSERC Discovery Grant to SB. 


REFERENCES 

Abramowitz M., Stegun I. A., 1964, Handbook of Math¬ 
ematical Functions with Formulas, Graphs, and Mathe¬ 
matical Tables, Dover, New York, NY 
Adams F. C., Fatuzzo M., 1996, ApJ, 464, 256 
Aitchison J., Brown J. A. C., 1957, The Lognormal Distri¬ 
bution, Cambridge University Press, London 
Allen A. P., Li B. L., Charnov E. L., 2001, Ecology Letters, 
4, 1 

Basu S., Jones C. E., 2004, MNRAS, 347, L47 
Bate M. R., Bonnell I. A., 2005, MNRAS, 356, 1201 
Chabrier G., 2003, PASP, 115, 763 

Chabrier G., 2005, in Corbelli E., Palla F., Zinnecker H., 
eds, Astrophysics and Space Science Library, Vol. 327, 
The Initial Mass Function 50 Years Later, Springer- 
Verlag, Dordrecht, p. 41 
Chabrier G., Baraffe I., 2000, ARA&A, 38, 337 
Clauset A., Shalizi C. R., Newman M. E. J., 2009, SIAM 
Review, 51, 661 

Crow E. L., Shimizu K., 1988, Lognormal Distributions: 

Theory and Applications, CRC Press, New York, NY 
Golberg M. A., 1984, An Introduction to Probability The¬ 
ory with Statistical Applications, Plenum Press, New 
York, NY 


© 0000 RAS, MNRAS 000, 000-000 























The MLP Distribution 7 


Gut A., 2005, Probability: A Graduate Course, Springer, 
New York, NY 

Hennebelle P., Chabrier G., 2009, ApJ, 702, 1428 
Kennicutt Jr R. C., 1998, in Gilmore, G., Howell, D., eds., 
ASP Conf. Ser. Vol. 142, The Stellar Initial Mass Func¬ 
tion, Astron. Soc. Pac., San Francisco, p. 1 
Kroupa P., 2001, MNRAS, 322, 231 
Kroupa P., 2002, Science, 295, 82 

Limpert E., Staliel W. A., Abbt M., 2001, BioScience, 51, 
341 

Mandelbrot B., 1997, Fractals and Scaling in Finance, 
Springer-Verlag, New York, NY 
Miller G. E., Scalo J. M., 1979, ApJS, 41, 513 
Mitzenmacher M., 2006, Internet Mathematics, 1, 226 
Myers P. C., 2009, ApJ, 706, 1341 
Myers P. C., 2014, ApJ, 781, 33 
Padoan P., Nordlund A., 2002, ApJ, 576, 870 
Pareto V., 1896, Cours d’conomie Politique, Droz, Geneva 
Reed W. J., 2002, Journal of Regional Science, 42, 1 
Reed W. J., 2003, Phys. A, 319, 469 
Salpeter E. E., 1955, ApJ, 121, 161 

Scalo J. M., 1998, in Gilmore, G., Howell, D., eds., ASP 
Conf. Ser. Vol. 142, The Stellar Initial Mass Function, 
Astron. Soc. Pac., San Francisco, p. 201 
Seggern D. H. V., 1990, CRC Handbook of Mathematical 
Curves and Surfaces, CRC Press, Boca Raton, FL 
Zipf G. K., 1949, Human Behavior and the Principle of 
Least Effort, Addison-Wesley, Reading, MA 


APPENDIX A: ERROR FUNCTION 

A1 Definition and Basic Properties 

The Error function, denoted by erf(x), and the Complemen¬ 
tary Error Function, denoted by erfc(x), are defined as 


erf(x) = —= 

\Ar 


/' 


e t dt , 


erfc(x) = 1 — erf(x) = —= 


r 

J X 


dt . 


(Al) 


(A2) 


Note from the definition that erf(x) is an odd function. Also, 
lim erf(x) = 1 (A3) 

x —>-oo 

(Abramowitz & Stegun 1964). The derivative and integral 
of erf(x) are: 


T x eri{x) = ji 6 "* 2 


(A4) 


/ 


erf(x)dx = xerf(a;) H — x + C , 


(A5) 


where C is an integration constant. Finally, the Taylor ex¬ 
pansion for erf(x) is given by 


erf(2; ) = -=£ 


(-ir* 


n 2n+l 


(2n + 1 )n\ Pit 


x 3 x 5 

x ~Y + w + - 


A plot of erfc(x) is presented in Figure Al. 


(A6) 



x 


Figure Al. Complementary error function, erfc(a;). 


A2 Related Integrals used in our analysis 

Consider the integral 

h = J exp[—(ax 2 + bx + c)]dx, a>0 (A7) 

By completing the square and letting u = yja [x + ^) 
it can be brought to 


h = ^ i/ — exp 
2 V a 


b 1 — 4ac 
4a 


erf ( pa 


x H- 

2a 


+ C , 
(A8) 


where C is an integration constant. In particular, using the 
properties of erf(x) and its relation to erfc(x), we have 


r° 2 

/ exp[— (ax + bx + c)]dx 
Jo 


1 /7r 


= 2 Va eXP 


b 2 — 4 ac 
4a 


erfc 


( 2 ^) • 


(A9) 


The result in Equation A8 can also be used to find in¬ 
tegrals of the form 


I 2 = J x exp [— (B + Clnx) 2 ] dx , 


(A10) 


where C ^ 0 and x > 0. By using the substitution v = 
B + C lnx, we can bring it to the form 


12 = c exp 


,2 . a aB 

~ {v +v c~-c- 


dv , 


which we can evaluate to 

' 0? + 4 01BC \ 


r pH 

12 = 2C exp 


4C 2 


erf ( B + C In x + 


2 CJ 


(All) 

+ C . 
(A12) 


APPENDIX B: DERIVATION OF THE 
RESULTS IN SECTION 5 

B1 The MLP as a Probability Density Function 

We can explicitly verify that the MLP density f(m) (Equa¬ 
tion 14) is indeed a valid probability density function for 
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all a > 0, i.e., it is positive and integrates to 1. Posi¬ 
tivity is clear from the definition. To verify the normal¬ 
ization condition, let A = (a/2) exp (apo + a 2 <To/2), B = 
(l/\/2) [aoo + (ito/oo)], and C = — l/(y/2ao). We then have 

p oo n oo 

/ f(m)dm = A / m l ' 1+a - 1 eric(B + Clnm)dm . (Bl) 
Jo Jo 

Letting dv = m~^ a+1 ^ dm and u = erfc(B + Clnm), we can 
write this integral as 



A 00 

- m “erfc(.B + Clnm) 

a 0 

exp [— ( B A Clnm) 2 ] dm . 


(B2) 


Consider the first term on the right hand side in Equa¬ 
tion B2. As m —> oo, erfc(B + Clnm) —> 2 (since C < 0 
and erf(a;) is odd), hence for a > 0 the term vanishes in this 
limit. As m —► 0, the limit takes on an indeterminate form. 
Applying L’Hospital’s rule we can see that the limit is also 
zero in this case: 


lim m “[erfc(B + Clnm)] 

m—>0 

2Q g — (B+C In m) 2 

=- -=— lim- 

y/lOtt ™-—>o m“ 

= —^ lim exp [— x 2 — ^-(x — B) 1 = 0 , (B3) 

L C J 

where x = B + Clnm. We can evaluate the second (integral) 
term with Equation A12: 


2CA 

A 
a 
2 A 


poo 

/ m -0+ 1 ) exp [— (B + Clnm) 2 ] dm 

Jo 


exp 


exp 


g? + 4 olBC \ 
~^C 2 ) 

a 2 + 4 olBC \ 

4 C 2 J 


x erf ^B + C In* + ^ j 


= 1 . 


o 

(B4) 


B2 Cumulative Distribution 


The cumulative distribution function, 

/ m 

f(t)dt , 

-OO 


(B5) 


and for our case /(m) is defined on [0, oo). To find the closed 
form we first we apply integration by parts, using also the 
same definitions for A, B, and C as in subsection Bl: 


F(m) 



t “erfc(B + Clnt) 


m 

0 


2AC 

aty/n 



t ~ (a+1 > 


exp [— ( B + Clnt) 2 ] dt . 


(B6) 


We again use the identity in Equation A12 to evaluate the 
integral term and obtain 


F(m) 


A 


a 


m “erfc(B + Clnm) 


A 

-exp 

a 



+ 4 aBC 
4 C 2 



B + C In m + 



(B7) 


which, upon returning to the original parameters, becomes 


F(m) = 2 er ^ c ( ~ 


ln(m) - [ip 
y/2ao 


- - exp ( afi o + 


„2 2 
a a 0 


"erfc 


/ aero 

VW 


ln(m) — no \ 
V2a 0 ) 
(B8) 


B3 Raw Moments 

Next we derive a closed form for arbitrary raw moments of 
the distribution, as well as an expression for its variance. 
Let M be an MLP random variable with probability density 
function f(m). The fcth raw moment of M, defined as the 
expectation value of M k , is given by 

POO 

E[M k ] = I m k f(m)dm 

Jo 

POO 

= A m fe ~i“+ 1 'erfc(B + Clnm)dm , (B9) 

Jo 

with A,B, and C as defined in subsection Bl. Before we 
arrive at a closed form expression for the moments we con¬ 
sider the convergence of the integral in Equation B9. This 
integral diverges for k ^ a. To see this, note that if k ^ a 
then p = a + l — k ^ 1. Now write Equation B9 as 


f 


eiicJB + Clnm) 


m v 


■L 

r 


a erfc (B + Clnm) 


m p 


dm 


erfc (B + Clnm) 


m p 


dm , 
(BIO) 


where a € (0, oo) such that erfc(R + Clnm) ^ 1, Vm ^ a. 
The existence of such a is ensured by the fact that for C < 0, 
erfc(B + Clnm) is a continuous strictly increasing positive 
function having an upper limit of 2. Then 


r°° 1 r°° 1 

/ -erfc(R + Clnm)dm ^ / - dm, (Bll) 

J a rn p J a m? 

where the integral on the right hand side diverges for p ^ 
1. Thus, the integral in Equation B9 diverges for k ^ a. 
Conversely, the moments are finite for k < a. For any a ^ 0, 


POO -J POO -j 

/ -erfc(B + Clnm)dm ^ 2 / - dm, (B12) 

J a m, p J a mV 

which converges for p = a — k + 1 > 1. Together with 
Equation BIO and Equation B3 this proves the existence 
of the moments. Suppose that a > k. Then by simply let¬ 
ting a —> a — fc>0in the integrand of Equation Bl, we 
can see from Equation B3 and Equation B4 that the fcth 
moment is given by 


. 9 4 

E\M k ] = exp 

a — k 


(a — k) 2 + 4(a — k)BC\ 


4C 2 


) 


a — k 


exp 


a 2 0 k 2 


+ Pok 


(B13) 
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B4 Mass Fraction 

We can obtain the expression for the mass fraction (up to a 
mass m), 

rm rm 

cf>(m ) = / m!f(m!)dm! = A / (m , ) _ “erfc(-B+Clnm , )dm , , 
Jo Jo 

using the same approaches employed in the computation of 
the cumulative distribution and the raw moments in the 
previous sections. Since the mean exists only if a > 1, it is 
also natural to consider the closed-form expression for the 
mass fraction for a > 1. This can be achieved by replacing 
all explicit appearances of a with a—1 in Equation B7, while 
maintaining the dependencies in the parameters 71(a), B(a) 
and C{a) in the original form, as we did for the calculation 
of the raw moments. Thus, 




(B14) 


which, upon returning to the original parameters, be¬ 
comes 



(B15) 


Note that, as expected, 4>{m) becomes the expression for the 
expectation value as m —> oo (for a > 1). 


© 0000 RAS, MNRAS 000, 000-000 












